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1.  Introduction 

Shock/Boundary  layer  interaction  is  still  an  important  problem  for  super- 
sonic aircraft  designers.  Such  phenomena  play  an  important  role  both  for 
internal  and  external  aerodynamic.  These  interactions  can  lead  to  an  in- 
crease of  drag,  separation,  loss  of  performances.  Moreover,  unsteadiness  of 
the  shock  produces  strong  constrains  on  the  structure,  noise  and  modifies 
local  heat  transfer. 

Since,  for  such  interactions,  results  obtained  with  the  RANS  approach 
are  not  satisfactory  [12],  it  appears  natural  to  assess  the  capacity  of  LES 
to  simulate  such  interaction.  From  a numerical  point  of  view,  it  necessi- 
tates the  use  of  schemes  designed  to  minimize  the  numerical  dissipation 
in  shock-free  region  of  the  flow  since  it  was  demonstrated  in  [1]  that  the 
numerical  dissipation  of  high-order  shock-capturing  scheme  can  exceed  the 
subgrid-scale  dissipation.  To  satisfy  this  requirement,  we  use  a strategy  [2] 
built  on  the  mixing  of  the  characteristic  based  filters  [3]  and  a sensor  [4] 
able  to  distinguish  a turbulent  fluctuation  from  a shock.  To  validate  this 
approach,  the  case  of  the  interaction  of  an  oblique  shock  with  a boundary 
layer  developing  on  a plane  plate  was  chosen  since  it  has  been  extensively 
studied  experimentally  at  IRPHE  [6]  [7]. 
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2.  Mathematical  model 

In  this  study,  we  have  used  the  system  of  filtered  Navier-Stokes  equations 
proposed  by  Vreman  [8].  Some  terms  are  neglected  as  in  [5].  The  subgrid- 
scale  model  is  the  selective  mixed  scale  model  proposed  by  Sagaut  [9].  This 
model  was  proven  to  behave  particularly  well  for  wall-bounded  flow. 

The  numerical  scheme  is  based  on  non-linear  WENO  filters  defined  in 
Ref.  [2] . Let  Un  denote  the  vector  of  the  conservative  variables  evaluated  at 
the  time  nAt  and  At  be  the  time  step,  [/(n+1)  the  vector  of  the  conservative 
variables  after  the  application  of  any  explicit  time  integration  scheme.  This 
vector  is  spatially  filtered  to  give  the  final  state  Un+1  (Un+l  = ^r(f7^n+1^)). 
The  main  point  is  that  the  time  advancement  is  performed  with  a non- 
dissipative  spatial  operator  (noted  L).  The  filtering  pass  is  decomposed  as 
follows: 

f/(n+i)  _ tf("+1))  = [Id  + pA tLf)(U^n+1})  (1) 

where  Lf  is  any  dissipative  operator,  Id  is  the  identity  and  the  switch  ft  is 
defined  as: 

/ P = 0 if  ® < 6 

\ 0 = 1 if*  >e  K ’ 

where  ft  is  the  sensor  of  Ducros  et  al.[ 4]  defined  as: 

* = (div(n))2 

(div(  u))2  + (rot(  u))2 

where  u denotes  the  velocity  vector.  This  sensor  was  demonstrated  to  be 
able  to  distinguish  a turbulent  fluctuation  from  a shock  in  Refs.  [4]  and  [2]. 

In  this  study,  the  time  integration  is  performed  by  means  of  a third-order 
accurate  TVD  Runge-Kutta  method  proposed  by  Shu  and  Osher  [10]:  Note 
that  L is  referred  to  as  ’’base  scheme”  and  can  be  any  qth.  order  accurate 
finite  volume  or  finite  difference  non-dissipative  scheme. 

As  mentioned  by  Yee  et  al.  [3],  Lf  can  be  the  dissipative  part  of  any 
shock-capturing  scheme.  Here,  such  operator  is  derived  of  ENO  schemes. 
The  dissipative  part  of  ENO  scheme  is  extracted  easily  by  subtracting  a 
centered  scheme  to  the  ENO  one. 

In  this  study,  we  have  used  a combination  of  a fourth-order  centered 
base  scheme  with  a fifth-order  accurate  WENO  filter.  The  threshold  e is 
fixed  to  0.04.  In  practice,  this  method  limits  the  computation  of  the  filter 
to  about  20  percents  of  the  total  number  of  grid  points. 

Viscous  fluxes  are  discretized  by  means  of  a second  order  accurate  cen- 
tered scheme.  The  CFL  number  is  fixed  to  0.5. 


LES  OF  SHOCK/BOUNDARY  LAYER  INTERACTION 


769 


Figure  1.  Computational  domain 


3.  Description  of  the  configuration 

The  shock  is  generated  by  a corner  fixed  to  the  upper  plate  of  the  wind 
tunnel  and  interacts  with  the  turbulent  boundary  layer  which  develops  on 
the  lower  plate  (which  is  assumed  to  be  adiabatic). 

The  length  of  the  computational  domain  was  restricted  to  the  measure- 
ment zone  (it  begins  at  x — 252  mm,  and  ends  at  x — 440  mm).  The 
height  of  the  domain  is  70,  7 mm  whereas  the  height  of  the  wind  tunnel  is 
120  mm.  The  computational  domain  is  presented  in  Fig.  1 (note  that  the 
axe  x denotes  the  longitudinal  direction,  y the  spanwise  direction  and  2 the 
wall-normal  direction).  The  width  of  the  computational  domain  is  about 
15  mm  in  the  spanwise  direction  (for  the  test  cases  A,  C and  D)  whereas 
the  wind  tunnel  is  10  times  larger. 

At  the  inflow  (at  x — 252  mm ),  the  temperature  Te  outside  of  the 
boundary  layer  is  144.6  K , the  density  is  equal  to  9.66  10~2  %/m3  and 
the  Mach  number  is  fixed  to  2.3.  The  Reynolds  number  Re  based  on  the 
displacement  length  <5i  = 3.535  mm  is  equal  to  19132.  Periodic  boundary 
conditions  are  applied  in  the  spanwise  direction  y.  Non  reflective  conditions 
are  applied  on  the  upper  frontier.  At  the  outflow,  a 13  mm  long  sponge 
zone  ensures  the  damping  of  turbulent  fluctuations  which  are  evacuated  by 
means  of  non-reflective  conditions.  No-slip  and  adiabaticity  conditions  are 
applied  on  the  wall. 

In  this  study,  four  simulations  have  been  carried  out  in  order  to  check 
the  sensitivity  of  the  results  to  computational  details.  The  characteristics  of 
these  simulations  are  reported  in  Table  1.  The  number  of  grid  points  (Nx, 
Ny , Nz)  in  each  direction  are  mentioned  and  the  size  of  the  computational 
meshes  (A+,  A+,  A+)  are  given  in  wall  units.  We  precise  that  these  wall 
units  are  based  on  a experimentally  measured  friction  velocity  of  24.75m/s 
at  x = 260  mm.  Test  case  A is  the  base  case.  In  test  case  B,  the  influence 
of  a doubling  of  the  domain  size  in  the  spanwise  direction  is  investigated. 
The  effect  of  a refined  grid  in  the  longitudinal  direction  is  studied  using 
case  C and  case  D is  devoted  to  the  investigation  of  the  importance  of  a 
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The  generation  of  realistic  inflow  conditions  is  a very  important  issue 
in  LES.  In  this  study,  we  follow  the  methodology  introduced  by  Lund  et 
al.  [11]  and  extended  to  compressible  flow  by  Urbin  et  al.  [13].  This  strategy 
prevents  the  thickening  of  the  boundary  layer  using  an  ad  hoc  re-scaling. 
The  loss  of  self-similarity  induced  by  the  shock  requires  the  Lund  procedure 
to  be  applied  on  a second  simulation  of  time-developing  canonical  turbulent 
boundary  layer.  A x — cste  plane  of  this  second  simulation  is  extracted  and 
introduced  in  the  shock/turbulence  LES  inflow  plane  at  each  time  step. 

4.  Analysis  of  the  results 

First,  iso- values  of  the  mean  pressure  are  presented  in  Fig.  2.  The  incident 
shock  (which  would  impact  near  x = 336  mm  in  absence  of  boundary  layer) 
curves  penetrating  the  boundary  layer  (thick  of  11  mm  at  the  inflow)  and 
reflects  on  the  sonic  line  (also  represented  on  the  picture)  as  an  expansion 
fan.  The  compression  related  to  the  rising  of  pressure  waves  in  the  subsonic 
part  of  the  boundary  layer  focuses  to  form  a reflected  shock  (whom  the 
trace  begins  in  the  vicinity  of  the  position  x = 290  mm  at  the  wall).  The 
continuation  of  this  shock  would  impact  on  the  wall  near  x = 275  mm  and 
marks  the  beginning  of  the  interaction  zone. 

Iso-values  of  longitudinal  velocity  fluctuations  are  shown  in  Fig.  3.  Up- 
stream the  interaction,  these  fluctuations  are  strong  in  the  near  wall  region. 
One  can  observe  that  the  fluctuations  are  amplified  by  a factor  two  under 
the  reflected  shock  (near  x = 290  mm).  The  explanation  given  by  Laurent 
[7]  to  explain  this  amplification  of  the  fluctuations  is  a linear  effect  by  rapid 
distorsion.  The  zone  concerned  by  a high  level  of  fluctuation  spreads  above 
the  separated  zone  (near  x = 320  mm).  In  agreement  with  experimental 
observations,  one  can  notice  an  alignment  of  the  isovalue  lines  of  the  fluc- 
tuations with  the  reflected  shock  just  downstream  of  this  one.  Downstream 
of  the  interaction  zone,  longitudinal  and  vertical  (not  shown)  fluctuations 
are  maximum  in  the  middle  of  the  boundary  layer  ( z/5  ~ 0,6).  However, 
one  can  observe  a second  extremum  of  longitudinal  fluctuations  in  the  near 
wall  region  close  to  the  outflow,  this  an  evidence  of  a (slow)  return  toward 
an  equilibrium  state.  A particular  attention  must  be  paid  to  the  interpre- 
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Figure  2.  Iso-values  of  the  mean  pressure  and  sonic  line  (in  black) 


Figure  3.  Iso-values  of  the  mean  longitudinal  velocity  fluctuations 


tation  of  data  concerning  the  fluctuations  in  the  interaction  zone  since  the 
unsteadiness  of  the  shocks  and  of  the  recirculation  bubble  can  be  at  the  ori- 
gin of  a part  of  the  fluctuations  (it  is  then  not  possible  to  consider  turbulent 
fluctuations). 

Longitudinal  evolution  of  the  displacement  thickness  is  displayed  in 
Fig.  4.  The  computation  of  this  integral  quantity  is  not  trivial  since  the 
external  velocity  Ue(x)  varies  in  the  longitudinal  direction.  The  potential 
flow  is  supposed  to  be  reached  when  the  Ducros  sensor  value  exceeds  0.1. 
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Figure  4-  Longitudinal  evolution  of  the  displacement  thickness 
Case  A , Case  B A,  Case  C , Case  D O,  HWA  O 


This  method  gives  results  coherent  with  measurement  since  the  value  of 
<5i  at  the  inflow  of  the  computational  domain  is  correctly  predicted  (at 
x = 260  mm)  after  a short  transient.  The  evolution  of  in  the  interaction 
zone  (290  mm  < x < 340  mm)  will  not  be  commented  in  absence  of  corre- 
sponding experimental  data.  Looking  at  figure  4,  one  can  observe  that  the 
displacement  thickness  is  multiplied  by  a ratio  larger  than  two  during  the 
interaction.  This  amplification  rate  is  well  predicted  by  computations  and 
during  the  relaxation  phase  where  <5i  decreases,  the  discrepancies  do  not 
exceed  10  %. 

The  longitudinal  evolution  of  the  friction  coefficient  Cf  is  presented  in 
Fig.  5 for  the  four  cases  and  compared  with  Hot  Wire  Anemometry  (HWA) 
data  [7].  At  the  inflow,  the  discrepancies  on  the  Cf  between  computations 
and  experience  is  lower  than  10  % for  cases  A and  B.  (such  kind  of  under- 
estimation is  classical  with  LES  [14]).  Cases  C (with  the  best  resolution) 
and  D give  levels  of  Cf  very  close  to  the  experiment.  The  flow  is  seen  to 
be  separated  for  cases  A,  B and  C between  x = 290  mm  and  x — 340  mm. 
The  quasi  totality  of  the  interaction  zone  is  concerned  by  the  separation. 
For  the  case  D (without  SGS  model),  the  velocity  fluctuations  in  the  near 
wall  zone  are  stronger  and  the  flow  separates  latter  (near  x = 305  mm). 
Just  after  the  interaction  zone  (x  = 340  mm),  experimental  evaluations 
based  on  hypothesis  of  the  existence  of  a logarithmic  zone  lead  to  a friction 
coefficient  clearly  positive  at  the  opposite  of  the  computations.  In  the  re- 
laxation zone,  the  increase  of  Cf  between  computations  and  experiment  is 
very  similar  and,  at  the  outflow,  the  discrepancies  are  comparable  to  those 
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Figure  5.  Longitudinal  evolution  of  friction  coefficient 

Case  A , Case  B A,  Case  C , Case  D O,  HWA  Q 


observed  at  the  inflow  (about  10  %).  The  results  quality  is  clearly  supe- 
rior to  the  one  obtained  with  RANS  computations  (see  [12])  in  the  same 
configuration  (but  at  Ma  = 2.9). 

Independently  of  the  analysis  of  the  results  concerning  physical  vari- 
ables, we  verify  in  Fig.  6 that  the  Ducros  sensor  applies  nearly  exclusively 
in  the  zone  where  we  have  noticed  in  Fig.  2 the  presence  of  shock  or  ex- 
pansion. Consequently,  the  SGS  model  is  effective  (his  effect  is  not  masked 
by  numerical  dissipation)  in  the  most  part  of  the  boundary  layer.  Never- 
theless, one  can  observe  low  values  of  the  sensor  on  the  upper  limit  of  the 
boundary  layer  where  the  adaptation  of  the  inflow  conditions  may  not  be 
perfect.  Even  if  this  strategy  of  minimization  of  the  numerical  dissipation 
would  merit  other  validation  (see  also  [2])  on  different  test  cases,  results 
presented  here  allows  to  be  confident  with  this  concept. 


5.  Conclusion 

Bidimensional  interaction  an  oblique  shock  with  a plane  plate  has  been 
studied  numerically  using  LES  and  compared  with  experimental  data.  Nu- 
merical results  are  in  quantitative  agreement  with  experimental  results  and 
LES  can  now  be  considered  as  a predictive  tool  for  such  physically  complex 
flow.  The  separated  zone  is  correctly  described.  The  effects  of  the  size  of 
the  domain  in  the  spanwise  direction,  of  the  resolution  in  the  longitudinal 
direction  and  the  presence  of  a SGS  model  do  not  appear  to  be  deciding. 
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Figure  6.  Iso-values  of  the  mean  PDF  of  Ducros  sensor 
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